home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_7.2 / fltrfreq / fltrfreq.c next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  13.2 KB  |  381 lines
  1. /* 
  2.  * fltrfreq.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* FLTRFREQ:    program performs lowpass, highpass, bandpass, or bandstop
  10.  *            filtering on image by applying filter in frequency domain
  11.  *                      usage: fltrfreq inimg outimg [-l f1] [-h f1] [-b f1 f2]
  12.  *                                      [-s f1 f2] [-n ORDER] [-g] [-p] [-I] [-L]
  13.  *        Note: Low/High-pass frequencies are specified as fraction of 1,
  14.  *              where 1 is normalized highest frequency.
  15.  *        Note: image sidelength must be power of 2; if not the case,
  16.  *              this program changes the size to be so, either the
  17.  *              next higher or lower power of 2 as user-chosen (<-p>).
  18.  *
  19.  */
  20.  
  21. #include <stdio.h>
  22. #include <math.h>
  23. #include <stdlib.h>
  24. #include <malloc.h>
  25. #include <string.h>
  26.  
  27. #include "ip.h"
  28.  
  29. #if defined(WIN32) || defined(LINUX)
  30. #define M_PI    3.14159265358979323846
  31. #define log2(a) (log(a)/ 0.6931471805599)
  32. #define exp2(a) (pow((double)2, (double)a))
  33. #endif
  34. #define PASS_TYPE_DFLT 0        /* default filter passband is low-pass */
  35. #define F1_DFLT 0.5             /* default filter cutoff is half */
  36. #define F2_DFLT 0.0
  37. #define ORDER_DFLT 1            /* default order of Butterworth filter */
  38. #define BUTTERWORTH 1           /* Butterworth filter type */
  39. #define GAUSSIAN 2              /* Gaussian filter type */
  40.  
  41. void input (int argc, char *argv[], short *passType, double *f1, double *f2, double *order, short *fltrType, short *zeroPad, long *invertFlag);
  42. void window (float **image, long nRow, long nCol, long nRowS, long nColS);
  43. void usage (short flag);
  44.  
  45. main (argc, argv)
  46.      int argc;
  47.      char *argv[];
  48.  
  49. {
  50.   register long x, y, nCol, nRow,  /* image input sidelengths */
  51.     nCol2, nRow2;               /* power of 2 spectrum sidelengths */
  52. //                  nColD2, nRowD2; /* middle row/column of spectrum image */
  53.   long xStart, yStart;          /* beginning coord.s of pow spec in image */
  54.   Image *imgI, *imgO;           /* I/O image structures */
  55.   unsigned char **imgIn,        /* image array */
  56.   **imgOut;                     /* output image */
  57.   short passType;               /* lowpass=0, high=1, band=2, stop=3 */
  58.   double f1, f2;                /* 3dB frequency bounds of filter */
  59.   double order;                 /* order of Butterworth filter */
  60.   short fltrType;               /* Butterworth or Gaussian filter type */
  61.   short zeroPad;                /* if 1, upsize image to pow of 2; 0 => down */
  62.   float **imgReal, **imgImag;   /* complex image arrays */
  63.   //float *wndwX, *wndwY;       /* row and column window vectors */
  64.   double min, max;              /* maximum value in power spectrum */
  65.   float max0, min0;             /* initial max,min values of image */
  66.   double norm;                  /* normalization factor for img intensities */
  67.   double temp;
  68.   long invertFlag;
  69.   //double  sqrt (), log2 (), log10 (), exp2 ();
  70.  
  71.   input (argc, argv, &passType, &f1, &f2, &order, &fltrType, &zeroPad, &invertFlag);
  72.  
  73. /* read input image */
  74.   imgI = ImageIn (argv[1]);
  75.   imgIn = imgI->img;
  76.   nRow = ImageGetHeight (imgI);
  77.   nCol = ImageGetWidth (imgI);
  78.   printf ("Original image is %3dx%3d.\n", nCol, nRow);
  79.  
  80. /* invert image intensities for going into processing */
  81.   if (invertFlag) {
  82.     for (y = 0; y < nRow; y++)
  83.       for (x = 0; x < nCol; x++)
  84.       imgIn[y][x] = 255 - imgIn[y][x];
  85.   }
  86.  
  87. /* check for image sidelengths powers of 2 */
  88.   nRow2 = nRow;
  89.   temp = log2 ((double) nRow);
  90.   if ((temp - (long) temp) != 0.0)
  91.     nRow2 = (zeroPad == 0) ?
  92.       (long) exp2 ((double) ((long) temp)) : (long) exp2 ((double) ((long) temp + 1));
  93.   nCol2 = nCol;
  94.   temp = log2 ((double) nCol);
  95.   if ((temp - (long) temp) != 0.0)
  96.     nCol2 = (zeroPad == 0) ?
  97.       (long) exp2 ((double) ((long) temp)) : (long) exp2 ((double) ((long) temp + 1));
  98.   if (nRow2 != nRow || nCol2 != nCol)
  99.     printf ("Image size changed to %dx%d so power of 2 (required for FFT).\n",
  100.             nRow2, nCol2);
  101.  
  102. /* allocate output image */
  103.   imgO = ImageAlloc (nRow2, nCol2, ImageGetDepth (imgI));
  104.   imgOut = ImageGetPtr (imgO);
  105.  
  106. /* allocate memory for complex array */
  107.   if ((imgReal = (float **) calloc (nRow2, sizeof (float *))) == NULL)
  108.       exit (2);
  109.   for (y = 0; y < nRow2; y++)
  110.     if ((imgReal[y] = (float *) calloc (nCol2, sizeof (float))) == NULL)
  111.         exit (3);
  112.  
  113.   if ((imgImag = (float **) calloc (nRow2, sizeof (float *))) == NULL)
  114.       exit (4);
  115.   for (y = 0; y < nRow2; y++)
  116.     if ((imgImag[y] = (float *) calloc (nCol2, sizeof (float))) == NULL)
  117.         exit (5);
  118.  
  119. /* center the image, and zero-pad or reduce size if required */
  120.   printf ("Forward 2-D FFT being performed.\n");
  121.   if (zeroPad == 0) {
  122.     yStart = (nRow - nRow2) / 2;
  123.     xStart = (nCol - nCol2) / 2;
  124.     for (y = 0, max0 = (float) 0.0, min0 = (float) 256.0; y < nRow2; y++) {
  125.       for (x = 0; x < nCol2; x++) {
  126.         imgReal[y][x] = (float) imgIn[y + yStart][x + xStart];
  127.         imgImag[y][x] = (float) 0.0;
  128.         if (imgReal[y][x] > max0)
  129.           max0 = imgReal[y][x];
  130.         if (imgReal[y][x] < min0)
  131.           min0 = imgReal[y][x];
  132.       }
  133.     }
  134.   }
  135.   else {
  136.     yStart = (nRow2 - nRow) / 2;
  137.     xStart = (nCol2 - nCol) / 2;
  138.     for (y = 0; y < nRow; y++) {
  139.       for (x = 0; x < nCol; x++) {
  140.         imgReal[y + yStart][x + xStart] = (float) imgIn[y][x];
  141.         imgImag[y][x] = (float) 0.0;
  142.       }
  143.     }
  144.   }
  145.  
  146. /* perform windowing */
  147.   if (zeroPad == 0)
  148.     window (imgReal, nRow2, nCol2, nRow2, nCol2);
  149.   else
  150.     window (imgReal, nRow2, nCol2, nRow, nCol);
  151.  
  152. /* perform 2-dimensional FFT */
  153.   fft2d (imgReal, imgImag, nRow2, nCol2, -1);
  154.  
  155. /* perform filtering on frequency-domain image */
  156.   printf ("Frequency-domain filtering being performed.\n");
  157.   if (fltrType == BUTTERWORTH)
  158.     fltrbttr (imgReal, imgImag, nRow2, nCol2, passType, f1, f2, order);
  159.   else
  160.     fltrgaus (imgReal, imgImag, nRow2, nCol2, passType, f1, f2);
  161.  
  162. /* perform inverse FFT */
  163.   printf ("Inverse 2-D FFT being performed.\n");
  164.   fft2d (imgReal, imgImag, nRow2, nCol2, 1);
  165.  
  166. /* find post-processing normalization factors for image */
  167.   for (y = 0, min = 1000.0, max = -1000.0; y < nRow2; y++) {
  168.     for (x = 0; x < nCol2; x++) {
  169.       if (imgReal[y][x] < min)
  170.         min = imgReal[y][x];
  171.       if (imgReal[y][x] > max)
  172.         max = imgReal[y][x];
  173.     }
  174.   }
  175.   if (min > 0.0)
  176.     min = 0.0;
  177.   norm = max0 / (max - min);
  178.  
  179.   for (y = 0; y < nRow2; y++)
  180.     for (x = 0; x < nCol2; x++)
  181.     imgOut[y][x] = (unsigned char) ((imgReal[y][x] - min) * norm + 0.5);
  182.  
  183. /* output image - un-invert it */
  184.   if (invertFlag) {
  185.     for (y = 0; y < nRow2; y++)
  186.       for (x = 0; x < nCol2; x++)
  187.       imgOut[y][x] = 255 - imgOut[y][x];
  188.   }
  189.  
  190. /* write output image */
  191.   ImageOut (argv[2], imgO);
  192.   return (0);
  193. }
  194.  
  195.  
  196. /* WINDOW:      function multiplies vector by smoothing window
  197.  *                    usage: window (image, nRow, nCol, nRowS, nColS)
  198.  *              Note: There are 2 sets of image sizes for the following
  199.  *              reason. The (nRow, nCol) size is the image size to be
  200.  *              windowed; that is, this is a power of 2 for the FFT.
  201.  *              The (nRowS, nColS) size is a smaller size that is the original
  202.  *              image size before zero-padding. This smaller image has
  203.  *              been placed in the middle of the zero-padded image. It is
  204.  *              the SMALLER IMAGE, NOT THE FINAL IMAGE that should be windowed.
  205.  */
  206.  
  207. void
  208. window (image, nRow, nCol, nRowS, nColS)
  209.      float **image;             /* image to be smoothed */
  210.      long nRow, nCol;           /* size of image to be windowed */
  211.      long nRowS, nColS;         /* size of smaller image b4 padding */
  212. {
  213.   float *hammingX, *hammingY;   /* horiz. and vert. hamming windows */
  214.   double nM1;                   /* no. rows/cols minus 1 */
  215.   long xStart, yStart;          /* offset smaller image in larger */
  216.   long x, y;
  217.  
  218. /* allocate space for hamming windows */
  219.   if ((hammingX = (float *) calloc (nCol, sizeof (*hammingX))) == NULL)
  220.     exit (1);
  221.   if ((hammingY = (float *) calloc (nRow, sizeof (*hammingY))) == NULL)
  222.     exit (2);
  223.  
  224. /* put window in middle of zeroed vector if image has been padded */
  225.   for (x = 0; x < nCol; x++)
  226.     hammingX[x] = (float) 0.0;
  227.   for (y = 0; y < nRow; y++)
  228.     hammingY[y] = (float) 0.0;
  229.   xStart = (nCol - nColS) / 2;
  230.   yStart = (nRow - nRowS) / 2;
  231.  
  232. /* calculate Hamming window */
  233.   nM1 = (double) (nColS - 1);
  234.   for (x = 0; x < nColS; x++)
  235.     hammingX[x + xStart] = (float) (0.54 - 0.46 * cos ((M_PI * 2.0 * (double) x) / nM1));
  236.   nM1 = (double) (nRowS - 1);
  237.   for (y = 0; y < nRowS; y++)
  238.     hammingY[y + yStart] = (float) (0.54 - 0.46 * cos ((M_PI * 2.0 * (double) y) / nM1));
  239.  
  240. /* apply window to image */
  241.   for (y = 0; y < nRow; y++)
  242.     for (x = 0; x < nCol; x++)
  243.       image[y][x] = image[y][x] * hammingX[x] * hammingY[y];
  244.  
  245. }
  246.  
  247.  
  248. /* USAGE:       function gives instructions on usage of program
  249.  *                      usage: usage (flag)
  250.  *              When flag is 1, the long message is given, 0 gives short.
  251.  */
  252.  
  253. void
  254. usage (flag)
  255.      short flag;                /* flag =1 for long message; =0 for short message */
  256. {
  257.  
  258. /* print short usage message or long */
  259.   printf ("USAGE: fltrfreq inimg outimg [-l f1] [-h f1] [-b f1 f2] \n");
  260.   printf ("                            [-s f1 f2] [-n ORDER] [-g] [-p] [-I] [-L]\n");
  261.  
  262.   if (flag == 0)
  263.     exit (1);
  264.  
  265.   printf ("\nfltrfreq performs frequency domain filtering on image.\n\n");
  266.   printf ("ARGUMENTS:\n");
  267.   printf ("    inimg: input image filename (TIF)\n");
  268.   printf ("   outimg: output image filename (TIF)\n\n");
  269.   printf ("OPTIONS:\n");
  270.   printf ("     -l f1: LOW-PASS FILTER passing frequencies below cutoff, f1.\n");
  271.   printf ("     -h f1: HIGH-PASS FILTER passing frequencies above cutoff, f1.\n");
  272.   printf ("  -b f1 f2: BAND-PASS FILTER passing frequencies between f1-f2.\n");
  273.   printf ("  -s f1 f2: STOP-BAND-PASS FILTER passing frequencies not f1-f2.\n");
  274.   printf ("  -n ORDER: order of Butterworth filter, dflt = %d.\n", ORDER_DFLT);
  275.   printf ("        -g: flag to perform Gaussian filter, default is Butterworth.\n");
  276.   printf ("        -p: flag for zero padding if set, and image row or column is\n");
  277.   printf ("            not a power of 2, image size is increased by zero-.\n");
  278.   printf ("            padding; otherwise image size is decreased (default).\n\n");
  279.   printf ("            NOTE -- The frequencies (f1,f2) should be expressed as a\n");
  280.   printf ("            number 0 to 1.0; this is a fractional frequency value of the\n");
  281.   printf ("            full original pass band. For example, a low-pass filter with\n");
  282.   printf ("            cutoff frequency of 0.5 will reduce the bandwidth by half.\n");
  283.   printf ("            Where the image is not square, the fractional frequency\n");
  284.   printf ("            value is relative to the higher frequency corresponding to\n");
  285.   printf ("            the longer x or y axis.\n");
  286.   printf ("        -I: invert image before processing\n");
  287.   printf ("        -L: print Software License for this module\n");
  288.  
  289.   exit (1);
  290. }
  291.  
  292.  
  293. /* INPUT:       function reads input parameters
  294.  *                   usage: input (argc, argv, &passType, &f1, &f2,
  295.  *                                      &order, &fltrType, &zeroPad)
  296.  */
  297.  
  298. void
  299. input (argc, argv, passType, f1, f2, order, fltrType, zeroPad, invertFlag)
  300.      int argc;
  301.      char *argv[];
  302.      short *passType;           /* lowpass=0, highpass=1, bandpass=2, stoppass= 3 */
  303.      double *f1, *f2;           /* 3db cutoff frequencies */
  304.      double *order;             /* order of Butterworth filter */
  305.      short *fltrType;           /* filter type */
  306.      short *zeroPad;            /* if 1, increase img size to 2-power; 0 => decrease */
  307.      long *invertFlag;          /* invert image before processing */
  308.  
  309. {
  310.   long n;
  311.   double atof ();
  312.  
  313.   if (argc < 3)
  314.     usage (1);
  315.  
  316.   *passType = PASS_TYPE_DFLT;
  317.   *f1 = F1_DFLT;
  318.   *f2 = F2_DFLT;
  319.   *order = ORDER_DFLT;
  320.   *fltrType = BUTTERWORTH;
  321.   *zeroPad = 0;
  322.   *invertFlag = 0;
  323.  
  324.   for (n = 3; n < argc; n++) {
  325.     if (strcmp (argv[n], "-l") == 0) {
  326.       if (++n == argc || argv[n][0] == '-')
  327.         usage (0);
  328.       *passType = 0;
  329.       *f1 = atof (argv[n]);
  330.     }
  331.     else if (strcmp (argv[n], "-h") == 0) {
  332.       if (++n == argc || argv[n][0] == '-')
  333.         usage (0);
  334.       *passType = 1;
  335.       *f1 = atof (argv[n]);
  336.     }
  337.     else if (strcmp (argv[n], "-b") == 0) {
  338.       if (++n == argc || argv[n][0] == '-')
  339.         usage (0);
  340.       *passType = 2;
  341.       *f1 = atof (argv[n]);
  342.       if (++n == argc || argv[n][0] == '-')
  343.         usage (0);
  344.       *f2 = atof (argv[n]);
  345.     }
  346.     else if (strcmp (argv[n], "-s") == 0) {
  347.       if (++n == argc || argv[n][0] == '-')
  348.         usage (0);
  349.       *passType = 3;
  350.       *f1 = atof (argv[n]);
  351.       if (++n == argc || argv[n][0] == '-')
  352.         usage (0);
  353.       *f2 = atof (argv[n]);
  354.     }
  355.     else if (strcmp (argv[n], "-n") == 0) {
  356.       if (++n == argc || argv[n][0] == '-')
  357.         usage (0);
  358.       *order = atof (argv[n]);
  359.     }
  360.     else if (strcmp (argv[n], "-g") == 0)
  361.       *fltrType = GAUSSIAN;
  362.     else if (strcmp (argv[n], "-p") == 0)
  363.       *zeroPad = 1;
  364.     else if (strcmp (argv[n], "-I") == 0)
  365.       *invertFlag = 1;
  366.     else if (strcmp (argv[n], "-L") == 0) {
  367.       print_sos_lic ();
  368.       exit (0);
  369.     }
  370.     else
  371.       usage (0);
  372.   }
  373.  
  374.   if (*f1 < 0.0 || *f1 >= 1.0 || *f2 < 0.0 || *f2 >= 1.0) {
  375.     printf ("Cutoff frequency values must be between 0.0 and 1.0.\n");
  376.     usage (1);
  377.   }
  378.  
  379.   return;
  380. }
  381.